# -*- eval: (visual-line-mode t); eval: (auto-fill-mode nil); fill-column: 1000000; org-ascii-text-width: 10000;  eval: (add-to-list 'org-odt-label-styles '("value0" "%c" "value0" "%n")); eval: (setq org-odt-category-map-alist '(("__Figure__" "Illustration" "value0" "Figure" org-odt--enumerable-image-p))); -*-
#+author: Николай Пузанов
#+email: punzik@gmail.com
#+language: ru
#+title: Быстрое вычисление медианы в целых числах
#+options: toc:nil

* Введение

#+begin_src elisp :results none :tangle no :exports none
  (if org-babel-restart-session
      (kill-system-buffer "Python"))
#+end_src

В задачах статистики и цифровой обработки сигналов нередко возникает потребность в вычислении медианы некоторого массива данных. Обычно нужно вычислить медиану от нескольких значений, что реализуется относительно просто и работает быстро. Но бывают случаи, когда массив содержит тысячи или десятки тысяч элементов, а машинного времени и ресурсов не очень много. В этом случае может помочь алгоритм binmedian со средней сложностью O(n), или его модификация binapprox, которая вычисляет приближенную медиану, но гарантированно за O(n).

В этой статье рассматривается рекурсивная реализация алгоритма binmedian, вычисляющая точное значение медианы для массива целых чисел конечной разрядности. Алгоритм имеет сложность строго O(n).

* Исходные данные

#+begin_src python :session :exports none :results none
  import random
  data_width = 5
  data_range = 2 ** data_width
  data_len = 15
#+end_src

В качестве примера возьмём массив из src_python[:session]{data_len} {{{results(=15=)}}} беззнаковых чисел (рис.1).

#+begin_src python :session :exports none :results table
  # numbers = [ random.randrange(0, data_range) for _ in range(data_len) ]
  numbers = [13, 2, 26, 31, 14, 10, 8, 28, 8, 18, 10, 13, 28, 31, 23]
  median_idx = round(len(numbers) / 2) - 1
  numbers
#+end_src

#+RESULTS:
| 13 | 2 | 26 | 31 | 14 | 10 | 8 | 28 | 8 | 18 | 10 | 13 | 28 | 31 | 23 |

#+caption: Рис.1: Исходный массив чисел. Верхняя строка - индексы элементов массива.
[[./pictures/p1.png]]

Найдём медиану обычным способом, т.е. отсортируем массив и возьмём средний элемент (рис.2).

#+begin_src python :session :exports none :results table
  numbers_sorted = numbers
  numbers_sorted.sort()
  numbers_sorted
#+end_src

#+RESULTS:
| 2 | 8 | 8 | 10 | 10 | 13 | 13 | 14 | 18 | 23 | 26 | 28 | 28 | 31 | 31 |

#+caption: Рис.2: Отсортированный массив. Средний элемент выделен пунктирным прямоугольником.
[[./pictures/p2.png]]


Средний элемент имеет номер src_python[:session]{median_idx} {{{results(=7=)}}}. Его значение, а соответственно и значение медианы равно src_python[:session]{numbers_sorted[median_idx]} {{{results(=14=)}}}. Это число мы будет использовать в дальнейшем для проверки корректности описанных методов.

* Метод binmedian
Как известно, сортировка - это довольно затратная задача для выполнения "в железе". К счастью, существует способ избежать этой процедуры с помощью метода вычисления медианы под названием binmedian. В общем случае он позволяет вычислить точное значение медианы без использования сортировки /в среднем/ за O(n). В случае целых чисел, как будет показано ниже, с помощью этого метода вычисляется точное значение /строго/ за O(n).

Для вычисления медианы методом binmedian нужно построить гистограмму входного массива чисел, а затем последовательно складывать значения бинов, пока сумма не станет больше половины длины последовательности (номера центрального элемента массива).

#+begin_src python :session :exports none :results table
  hist_len = 4

  nmin = min(numbers)
  nmax = max(numbers)
  hist_bins = [nmin] + [round(x * ((nmax-nmin)/hist_len)) for x in range(1, hist_len)] + [nmax]

  hist_bins = []
  bin_gap = (nmax - nmin + 1) / hist_len

  for n in range(0, hist_len):
      ndn = round(n * bin_gap + nmin)
      nup = round((n + 1) * bin_gap + nmin)
      hist_bins.append((ndn, nup))

  hist_bins_str = ", ".join(["[{};{})".format(b[0], b[1]) for b in hist_bins])

  [["Номер бина:"] + list(range(len(hist_bins))),
   ["Интервал:"] + ["[{};{})".format(b[0], b[1]) for b in hist_bins]]
#+end_src

#+RESULTS:
| Номер бина: |      0 |       1 |       2 |       3 |
| Интервал:   | [2;10) | [10;17) | [17;24) | [24;32) |

В качестве примера построим гистограмму из четырёх бинов. Для этого разобьем диапазон чисел в массиве на четыре части. С учётом входных значений, диапазоны бинов будут следующие (каждое значение - это полуинтервал [a;b), в него входят значения большие или равные левой границе, и меньшие правой границы): src_python[:session]{hist_bins_str} {{{results(=[2;10)\, [10;17)\, [17;24)\, [24;32)=)}}}.

На рис.3 цветами выделены интервалы в отсортированном массиве.

#+caption: Рис.3: Отсортированный массив, разбитый на интервалы.
[[./pictures/p3.png]]

#+begin_src python :session :exports none :results table
  hist = [0] * hist_len

  for x in numbers:
      for n, b in enumerate(hist_bins):
          if (x >= b[0] and x < b[1]):
              hist[n] += 1

  [["Номер бина:"] + list(range(len(hist))),
   ["Кол-во вхождений:"] + hist]
#+end_src

#+RESULTS:
| Номер бина:       | 0 | 1 | 2 | 3 |
| Кол-во вхождений: | 3 | 5 | 2 | 5 |

#+begin_src python :session :exports none :results value
  bin_sum = 0
  bin_median_idx = -1

  for n, h in enumerate(hist):
      bin_sum += h
      if (bin_sum > median_idx):
          bin_median_idx = n
          break

  bin_median = hist_bins[bin_median_idx]
  bin_median_str = "[{};{})".format(bin_median[0], bin_median[1])
  bin_median_str
#+end_src

#+RESULTS:
: [10;17)

Вычислим значения бинов гистограммы, а затем просуммируем их слева направо до тех пор, пока сумма не превысит число src_python[:session]{median_idx} {{{results(=7=)}}} (индекс среднего элемента массива). Бин, на котором остановился счёт будет тем бином, в интервале которого находится медиана. В нашем случае это бин под номером src_python[:session]{bin_median_idx} {{{results(=1=)}}}, т.е. медиана находится в интервале src_python[:session]{bin_median_str} {{{results(=[10;17)=)}}}. Гистограмма показана на рис.4.

#+caption: Рис.4. Гистограмма. В нижней части обозначены номера бинов, затем выше интервалы бинов, и наконец их значения. Сумма первых двух бинов равна src_python[:session]{bin_sum} {{{results(=8=)}}}, что больше, чем индекс среднего элемента массива.
[[./pictures/p4.png]]

Нетрудно вычислить сложность этого метода. Количество операций равно n * k + c, где n - количество элементов в массиве, k - количество операций, необходимых для обновления гистограммы, c - расходы на подготовку гистограммы и сложение значений бинов. Т.к. k и c константы (c зависит от разрядности чисел и в общем случае не зависит от количества элементов), то время вычисления зависит (линейно) только от количества элементов массива. Таким образом сложность метода равняется O(n).

* Точное вычисление медианы
Чтобы вычислить точное значение медианы, очевидно нужно сделать гистограмму такой длины, чтобы интервал каждого бина был равен единице. Тогда результирующее значение точно попадёт в значение медианы (напомню, мы производим вычисления в целых числах).

В нашем случае необходима гистограмма от минимального значения (src_python[:session]{nmin} {{{results(=2=)}}}) до максимального (src_python[:session]{nmax} {{{results(=31=)}}}) с количеством бинов src_python[:session]{nmax-nmin+1} {{{results(=30=)}}}.

Чтобы упростить схему (а мы же в итоге планируем сделать это в железе?), можно взять полный интервал разрядности чисел. В этом случае адресом бина будет служить просто значение очередного числа из массива. Для нашего примера я специально выбрал числа шириной src_python[:session]{data_width} {{{results(=5=)}}} бит, то есть их значения лежат в диапазоне от =0= до src_python[:session]{data_range-1} {{{results(=31=)}}}.

#+begin_src python :session :exports none :results table
    hist_full = [0] * data_range

    for n in numbers:
        hist_full[n] += 1

    [["Номер бина:"] + list(range(len(hist_full))),
     ["Кол-во вхождений:"] + hist_full]
    # [["Номер бина", "Кол-во вхождений:"]] + list(enumerate(hist_full))
#+end_src

#+RESULTS:
| Номер бина:       | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 |
| Кол-во вхождений: | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 2 | 0 |  2 |  0 |  0 |  2 |  1 |  0 |  0 |  0 |  1 |  0 |  0 |  0 |  0 |  1 |  0 |  0 |  1 |  0 |  2 |  0 |  0 |  2 |

#+begin_src python :session :exports none :results none
  bin_sum = 0
  median = -1

  for n, h in enumerate(hist_full):
      bin_sum += h
      if (bin_sum > median_idx):
          median = n
          break
#+end_src

Построим гистограмму. Как и в предыдущем случае, будем складывать значения бинов слева направо пока сумма не превысит src_python[:session]{median_idx} {{{results(=7=)}}} (рис.5). В итоге мы остановимся на бине src_python[:session]{median} {{{results(=14=)}}} с суммой, равной src_python[:session]{bin_sum} {{{results(=8=)}}}. Номер бина и есть значение медианы.

#+caption: Рис.5. Гистограмма с количеством бинов, равном полному диапазону входных чисел. Индекс выделенного бина и есть искомая медиана.
[[./pictures/p5.png]]

Сложность алгоритма не изменилась, но увеличился размер гистограммы. А что если у нас будут числа большей разрядности? Строить гистограмму например для 32-битных данных очень расточительное занятие - понадобится массив памяти на 4 миллиарда слов.

* Точное вычисление медианы для чисел с большой разрядностью
В этом случае можно разбить вычисление на несколько этапов. Сначала, используя старшие разряды чисел вычислить приближенное значение медианы. Затем, опускаясь по разрядам, прийти к точному значению. Количество этапов может быть произвольным в зависимости от разрядности чисел и требований к объёму используемой памяти.

#+begin_src python :session :exports none :results none
  high_width = round(data_width / 2 + 0.5)
  low_width = data_width - high_width
#+end_src

У нас разрядность небольшая, по этому вычислим медиану в два этапа. Для этого разобьём наш интервал разрядности (src_python[:session]{data_width} {{{results(=5=)}}} бит) на две части. Первая часть будет состоять из src_python[:session]{high_width} {{{results(=3=)}}} бит, вторая - из src_python[:session]{low_width} {{{results(=2=)}}} бит.

Построим гистограмму длиной src_python[:session]{2 ** high_width} {{{results(=8=)}}} используя только src_python[:session]{high_width} {{{results(=3=)}}} старших бит для индексации бинов.

#+begin_src python :session :exports none :results table
  high_hist = [0] * (2 ** high_width)

  for x in numbers:
      idx = x >> low_width
      high_hist[idx] += 1

  [["Номер бина:"] + list(range(len(high_hist))),
   ["Кол-во вхождений:"] + high_hist]
#+end_src

#+RESULTS:
| Номер бина:       | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| Кол-во вхождений: | 1 | 0 | 4 | 3 | 1 | 1 | 1 | 4 |

#+begin_src python :session :exports none :results none
  bin_sum = 0
  median_high = -1

  for n, h in enumerate(high_hist):
      bin_sum += h
      if (bin_sum > median_idx):
          median_high = n
          break
#+end_src

Так же, как и в предыдущих случаях, просуммируем бины слева направо и определим номер бина, на котором сумма превысит индекс центрального элемента. В нашем случае это бин с номером src_python[:session]{median_high} {{{results(=3=)}}} (рис.6). Значение старших src_python[:session]{high_width} {{{results(=3=)}}} бит медианы будет равно этому номеру, что в бинарном виде будет выглядеть как src_python[:session]{"0b{0:0{len}b}".format(median_high, len = high_width)} {{{results(=0b011=)}}}.

#+caption: Рис.6. Гистограмма, посроенная по старшим src_python[:session]{high_width} {{{results(=3=)}}} битам входных чисел. Соответствующие старшие биты медианы равны индексу бина, на котором остановился подсчёт суммы.
[[./pictures/p6.png]]

Теперь можно вычислить точное значение медианы, используя наше знание о значении её старших бит. Для этого построим вторую гистограмму. Для её построения будем отбирать числа, старшие биты которых равны значению, полученному на предыдущем шаге - src_python[:session]{"0b{0:0{len}b}".format(median_high, len = high_width)} {{{results(=0b011=)}}}. Для индексации гистограммы будем использовать только младшие src_python[:session]{low_width} {{{results(=2=)}}} бит этих чисел, т.е. длина гистограммы будет равняться src_python[:session]{2 ** low_width} {{{results(=4=)}}}.

На этом шаге мы должны ввести новую операцию, которой не было в ранее. Необходимо запомнить количество чисел, значение старших бит которых меньше значения, полученного на предыдущем шаге. Фактически это ещё один бин гистограммы с интервалом "меньше меньшего". Значение этого бина будет необходимо на последнем шаге вычислений.

#+begin_src python :session :exports none :results table
  low_hist = [0] * (2 ** low_width)
  mask = (low_width ** 2) - 1
  left_sum = 0

  for x in numbers:
      h_bits = x >> low_width

      if (h_bits < median_high):
          left_sum += 1
      elif (h_bits == median_high):
          idx = x & mask
          low_hist[idx] += 1

  [["Номер бина:"] + list(range(len(low_hist))),
   ["Кол-во вхождений:"] + low_hist]
#+end_src

#+RESULTS:
| Номер бина:       | 0 | 1 | 2 | 3 |
| Кол-во вхождений: | 0 | 2 | 1 | 0 |

Остался последний шаг - просуммировать бины и остановиться на том, где сумма превысит индекс центрального элемента. Однако, в отличие от предыдущих примеров, начальное значение суммы равно не нулю, а количеству чисел, старшие биты которых меньше старших битов медианы (тот самый бин "меньше меньшего"). Нетрудно посчитать их количество - src_python[:session]{left_sum} {{{results(=5=)}}} штук.

#+begin_src python :session :exports none :results none
  bin_sum = left_sum
  median_low = -1

  for n, h in enumerate(low_hist):
      bin_sum += h
      if (bin_sum > median_idx):
          median_low = n
          break
#+end_src

#+begin_src python :session :exports none :results none
  median = (median_high << low_width) | median_low
#+end_src

В результате последовательного суммирования (показано на рис.7) мы получили номер бина - src_python[:session]{median_low} {{{results(=2=)}}}. Это и есть значение младших бит искомой медианы - src_python[:session]{"0b{0:0{len}b}".format(median_low, len = low_width)} {{{results(=0b10=)}}}. Склеив всё вместе получим ответ - src_python[:session]{"0b{0:0{len}b}".format(median, len = data_width)} {{{results(=0b01110=)}}}, или src_python[:session]{median} {{{results(=14=)}}}.

#+caption: Рис.7. Гистограмма, посроенная по младшим src_python[:session]{low_width} {{{results(=2=)}}} битам входных чисел. На рисунке выделен бин, индекс которого равен младшим битам значения медианы.
[[./pictures/p7.png]]

Несмотря на то, что число итераций удвоилось, время потраченное на вычисления всё ещё линейно зависит от количества элементов массива. Т.е. сложность алгоритма осталась той же.

* Заключение
В статье был показан метод быстрого вычисления медианы для массива челых чисел с использованием алгоритма binmedian. Автором был реализован описанный алгоритм в виде RTL, который был успешно использован в ПЛИС Xilinx серий 7 и US+.

Исходный код статьи вы сможете найти по ссылке: https://github.com/punzik/fast-median-article. Там вы так же найдёте код на Python, которым вычислялись значения, приведенные в статье.
